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Abstract 

I review several issues related to statistical description of gravitating 
systems in both static and expanding backgrounds. After briefly review- 
ing the results for the static background, I concentrate on gravitational 
clustering of collision-less particles in an expanding universe. In particu- 
lar, I describe (a) how the non linear mode-mode coupling transfers power 
from one scale to another in the Fourier space if the initial power spec- 
trum is sharply peaked at a given scale and (b) what are the asymptotic 
characteristics of gravitational clustering that are independent of the ini- 
tial conditions. Numerical simulations as well as analytic work shows 
that power transfer leads to a universal power spectrum at late times, 
somewhat reminiscent of the existence of Kolmogorov spectrum in fluid 
turbulence. 



1 Overview of the key issues and results 

The statistical mechanics of systems dominated by gravity has close connec- 
tions with areas of condensed matter physics, fluid mechanics, re-normalization 
group etc. and poses an incredible challenge as regards the basic formulation. 
The ideas also find application in many different areas of astrophysics and cos- 
mology, especially in the study of globular clusters, galaxies and gravitational 
clustering in the expanding universe. (For a overall review of statistical me- 
chanics of gravitating systems, see [I], [2], [3]; review of gravitational clustering 
in expanding background is available in [3] and in several textbooks in cosmol- 
ogy [5]; for a sample of different attempts to understand these phenomena by 
different groups; see [S], [7], [S], [5], [TU], [H] and the references cited therein.) 
It will be useful to begin with a broad overview and a description of the issues 
which will be addressed in this article. 

In Newtonian theory, the gravitational force can be described as a gradient 
of a scalar potential and the evolution of a set of particles under the action of 
gravitational forces can be described the equations 

K^ = -\7(f>{M,t); V2(/)=-47rG^m,(5c(x-x,) (1) 
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where is the position of the i—th particle, rrii is its mass. For an isolated 
system with sufBcicntly large number of particles, it is useful to investigate 
whether some kind of statistical description of such a system is possible. Such 
a description, however, is complicated by the long range, unscreened, nature of 
gravitational force. If a self gravitating system is divided into two parts, the 
total energy of the system cannot be expressed as the sum of the gravitational 
energy of the components. The conventional results in statistical physics are 
based on the extensivity of the energy which is clearly invalid for gravitating 
systems. To construct the statistical description of such a system, one must 
begin with the construction of the micro-canonical ensemble describing such a 
system. If the Hamiltonian of the system is H{pi,qi) then the volume g{E) of 
the constant energy surface H{pi,qi) — E will be of primary importance in the 
micro-canonical ensemble. The logarithm of this function will give the entropy 
S{E) = \ng{E) and the temperature of the system will be T{E) = (3{E)-^ = 
{dS/dE)-\ 

Systems for which a description based on canonical ensemble is possible, the 
Laplace transform of g{E) with respect to a variable (3 will give the partition 
function Z{[3). It is, however, trivial to show that gravitating systems of interest 
in astrophysics cannot be described by a canonical ensemble [T] , [T^] , [H] . Virial 
theorem holds for such systems and we have {2K+U) = where K and U are the 
total kinetic and potential energies of the system. This leads to E = K + U = 
~K; since the temperature of the system is proportional to the total kinetic 
energy, the specific heat will be negative: Cy = {dE/dT)v oc {dE/dK) < 0. 
On the other hand, the specific heat of any system described by a canonical 
ensemble Cy — 0^ < (A£')^ > will be positive definite. Thus one cannot 
describe self gravitating systems of the kind we are interested in by canonical 
ensemble. 

One may attempt to find the equilibrium configuration for self gravitating 
systems by maximizing the entropy S{E) or the phase volume g{E). It is again 
easy to show that no global maximum for the entropy exists for classical point 
particles interacting via Newtonian gravity. To prove this, we only need to 
construct a configuration with arbitrarily high entropy which can be achieved 
as follows: Consider a system of N particles initially occupying a region of finite 
volume in phase space and total energy E. We now move a small number of 
these particles (in fact, a pair of them, say, particles 1 and 2 will do) arbitrarily 
close to each other. The potential energy of interaction of these two particles, 
—Gmim2/ri2, will become arbitrarily high as ri2 — » 0. Transferring some of 
this energy to the rest of the particles, we can increase their kinetic energy 
without limit. This will clearly increase the phase volume occupied by the 
system without bound. This argument can be made more formal by dividing 
the original system into a small, compact core and a large diffuse halo and 
allowing the core to collapse and transfer the energy to the halo. 

The absence of the global maximum for entropy — as argued above — de- 
pends on the idealization that there is no short distance cut-off in the interaction 
of the particles, so that we could take the limit ri2 — *■ 0. If we assume, instead, 
that each particle has a minimum radius a, then the typical lower bound to 
the gravitational potential energy contributed by a pair of particles will be 
—Gmim2/2a. This will put an upper bound on the amount of energy that can 
be made available to the rest of the system. 

We have also assumed that part of the system can expand without limit — in 
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the sense that any particle with sufficiently large energy can move to arbitrarily 
large distances. In real life, no system is completely isolated and eventually one 
has to assume that the meandering particle is better treated as a member of 
another system. One way of obtaining a truly isolated system is to confine the 
system inside a spherical region of radius R with, say, reflecting wall. (Most 
of our discussion is confined to 3-dimensions and the situation is diffrent in 
2-dimensions; see e.g [TJ) 

The two cut-offs a and R will make the upper bound on the entropy finite, but 
even with the two cut-offs, the primary nature of gravitational instability cannot 
be avoided. The basic phenomenon described above (namely, the formation of 
a compact core and a diffuse halo) will still occur since this is the direction of 
increasing entropy. Particles in the hot diffuse component will permeate the 
entire spherical cavity, bouncing off the walls and having a kinetic energy which 
is significantly larger than the potential energy. The compact core will exist as a 
gravitationally bound system with very little kinetic energy. A more formal way 
of understanding this phenomena is based on the virial theorem for a system 
with a short distance cut-off confined to a sphere of volume V. In this case, the 
virial theorem will read as [2] 

2T+U = 3PV + $ (2) 

where P is the pressure on the walls and $ is the correction to the potential 
energy arising from the short distance cut-off. This equation can be satisfied 
in essentially three different ways. If T and U are significantly higher than 
3PV, then we have 2T + U ^ which describes a self gravitating systems in 
standard virial equilibrium but not in the state of maximum entropy. If T S> ?7 
and 3PV ^ one can have 2T w 3PV which describes an ideal gas with no 
potential energy confined to a container of volume V; this will describe the hot 
diffuse component at late times. If T <C [/ and 3PV <C $, then one can have 
t/ ~ $ describing the compact potential energy dominated core at late times. 
In general, the evolution of the system will lead to the production of the core 
and the halo and each component will satisfy the virial theorem in the form 
Such an asymptotic state with two distinct phases [IS] is quite different from 
what would have been expected for systems with only short range interaction. 
Considering its importance, I shall briefly describe in section [5] a toy model 
which captures the essential physics of the above system in an exactly solvable 
context. 

The above discussion focussed on the existence of global maximum to the 
entropy and we proved that it does not exist in the absence of two cut-offs. 
It is, however, possible to have local extrema of entropy which are not global 
maxima. Intuitively, one would have expected the distribution of matter in the 
configuration which is a local extrema of entropy to be described by a Boltzmann 
distribution, with the density given by p{x.) (x exp[— /?0(x)] where is the 
gravitational potential related to p by Poisson equation. This is indeed true; for 
a formal proof see [1] . This configuration is usually called the isothermal sphere 
(because it can be shown that, among all solutions to this equation, the one 
with spherical symmetry maximizes the entropy) and it is a local maximum of 
entropy. The second (functional) derivative of the entropy with respect to the 
configuration variables will determine whether the local extremum of entropy is 
a local maximum or a saddle point [16j , [17j . 



3 



The relevance of the long range of gravity in all the above phenomena can 
be understood by studying model systems with an attractive potential vary- 
ing as r~" with different values for a. Such studies confirm the results and 
interpretation given above; (see [IHI and references cited therein). 

Let us now consider the situation in the context of an expanding background. 
There is considerable amount of observational evidence to suggest that one of 
the dominant energy densities in the universe is contributed by self gravitating 
point particles. The smooth average energy density of these particles drive 
the expansion of the universe while any small deviation from the homogeneous 
energy density will cluster gravitationally. [For a review of cosmology from 
a contemperorary perspective, see e.g., [19] ] One of the central problems in 
cosmology is to describe the non linear phases of this gravitational clustering 
starting from a initial spectrum of density fluctuations. It is often enough 
(and necessary) to use a statistical description and relate different statistical 
indicators (like the power spectra, nth order correlation functions ....) of the 
resulting density distribution to the statistical parameters (usually the power 
spectrum) of the initial distribution. The relevant scales at which gravitational 
clustering is non linear are less than about 10 Mpc (where 1 Mpc = 3 x 10^'* cm 
is the typical separation between galaxies in the universe) while the expansion 
of the universe has a characteristic scale of about few thousand Mpc. Hence, 
non linear gravitational clustering in an expanding universe can be adequately 
described by Newtonian gravity provided the rescaling of lengths due to the 
background expansion is taken into account. This is easily done by introducing 
a proper coordinate for the j— th particle r^, related to the comoving coordinate 
Xi, by Yi — a{t)xi with a{t) describing the stretching of length scales due to 
cosmic expansion. The Newtonian dynamics works with the proper coordinates 
Ti which can be translated to the behaviour of the comoving coordinate 
by this rescaling. [This implies that, for all practical purposes, we are still in 
the domain of Newtonian gravity. There is a far deeper connection between 
thermodynamics and gravity [20] in the general relativistic domain which we 
will not discuss in these lectures.] 

As to be expected, cosmological expansion completely changes the nature 
of the problem because of several new factors which come in: (a) The problem 
has now become time dependent and it will be pointless to look for equilibrium 
solutions in the conventional sense of the word, (b) On the other hand, the 
expansion of the universe has a civilizing influence on the particles and acts 
counter to the tendency of gravity to make systems unstable, (c) In any small 
local region of the universe, one would assume that the conclusions describing a 
finite gravitating system will still hold true approximately. In that case, particles 
in any small sub region will be driven towards configurations of local extrema 
of entropy (say, isothermal spheres) and towards global maxima of entropy (say, 
core-halo configurations). 

An extra feature comes into play as regards the expanding halo from any sub 
region. The expansion of the universe acts as a damping term in the equations 
of motion and drains the particles of their kinetic energy — which is essentially 
the lowering of temperature of any system participating in cosmic expansion. 
This, in turn, helps gravitational clustering since the potential wells of nearby 
sub regions can capture particles in the expanding halo of one region when the 
kinetic energy of the expanding halo has been sufficiently reduced. 

The actual behaviour of the system will, of course, depend on the form of 
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a(t). However, for understanding the nature of clustering, one can take a{t) cx 
^2/3 -^Jiicii describes a matter dominated universe with critical density. Such 
a power law has the advantage that there is no intrinsic scale in the problem. 
Since Newtonian gravitational force is also scale free, one would expect some 
scaling relations to exist in the pattern of gravitational clustering. Converting 
this intuitive idea into a concrete mathematical statement turns out to be non 
trivial and difficult. 

It is clear that cosmological expansion introduces several new factors into 
the problem when compared with the study of statistical mechanics of isolated 
gravitating systems. (For a general review of statistical mechanics of gravitating 
systems, see [1]. For a sample of different approaches, see [21] and the references 
cited therein. Review of gravitational clustering in expanding background is also 
available in several textbooks in cosmology Though this problem can 

be tackled in a 'practical' manner using high resolution numerical simulations 
(for a review, see [23j). such an approach hides the physical principles which 
govern the behaviour of the system. To understand the physics, it is necessary 
to attack the problem from several directions using analytic and semi analytic 
methods. Several such attempts exist in the literature based on Zeldovich(like) 
approximations [23], path integral and perturbative techniques [SF], nonlinear 
scaling relations [26j and many others. In spite of all these it is probably fair 
to say that we still do not have a clear analytic grasp of this problem, mainly 
because each of these approximations have different domains of validity and do 
not arise from a central paradigm. 

I propose to attack the problem from a different angle, which has not re- 
ceived much attention in the past. The approach begins from the dynamical 
equation for the the density contrast in the Fourier space and casts it as an 
integro-differential equation. Though this equation is known in the literature 
(see, e.g. [12 )j received very little attention because it is not 'closed' 

mathematically; that is, it involves variables which are not natural to the for- 
malism and thus further progress is difficult. I will, however, argue that there 
exists a natural closure condition for this equation based on Zeldovich approxi- 
mation thereby allowing us to write down a closed integro-differential equation 
for the gravitational potential in the Fourier space. 

It turns out that this equation can form the basis for several further inves- 
tigations some of which are described in ref. |27| and in the second reference in 
[T]. Here I will concentrate on just two specific features, centered around the 
following issues: 

• If the initial power spectrum is sharply peaked in a narrow band of wave- 
lengths, how does the evolution transfer the power to other scales? In 
particular, does the non linear evolution in the case of gravitational inter- 
actions lead to a universal power spectrum (like the Kolmogorov spectrum 
in fluid turbulence)? 

• What is the nature of the time evolution at late stages? Does the gravita- 
tional clustering at late stages wipe out the memory of initial conditions 
and evolve in a universal manner? 

Fair amount of progress can be made as regards these questions using the 
integro-differential equation mentioned above and some of these aspects will 
be discussed in detail. 
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2 Phases of the self gravitating system 



As described in section [T] the statistical mechanics of finite, self gravitating, 
systems have the following characteristic features: (a) They exhibit negative 
specific heat while in virial equilibrium, (b) They are inherently unstable to 
the formation of a core-halo structure and global maximum for entropy does 
not exist without cut-offs at short and large distances, (c) They can be broadly 
characterized by two phases — one of which is compact and dominated by 
potential energy while the other is diffuse and behaves more or less like an ideal 
gas. The purpose of this section is to describe a simple toy model which exhibits 
all these features and mimics a self gravitating system [1] . 

Consider a system with two particles described by a Hamiltonian of the form 

77(P,Q;p,r) = ^ + H!„^ (3) 
2M 2/i r 

where (Q, P) are coordinates and momenta of the center of mass, (r, p) are the 
relative coordinates and momenta, M — 2m is the total mass, 11 — m/2 is the 
reduced mass and m is the mass of the individual particles. This system may 
be thought of as consisting of two particles (each of mass m) interacting via 
gravity. We shall assume that the quantity r varies in the interval {a,R). This 
is equivalent to assuming that the particles are hard spheres of radius a/2 and 
that the system is confined to a spherical box of radius R. We will study the 
"statistical mechanics" of this simple toy model. 

To do this, we shall start with the volume g{E) of the constant energy surface 
H = E. Straightforward calculation gives 



g{E) = AR" 





r^dr 




r 


E+ 


J a 




r 



(4) 



where A = (647r^m'^/3). The range of integration in ^ should be limited to the 
region in which the expression in the square brackets is positive. So we should 
use rmax = if {-Gm^/a) < E < {-Gm'^/R), and use r^ax = i? if 

{—Gm^/R) < E < +00. Since H > {—Gm? /a), we trivially have g{E) = for 
E < {—Gm?/a). The constant A is unimportant for our discussions and hence 
will be omitted from the formulas hereafter. The integration in (jj]) gives the 
following result: 



9{E) 
[Gm^f 



^i-E)-^ (1 + , i-Gmya) <E< {-Gm^ / R) 

RE\^ , aE n3 



-i-E)-^ [(1 + - (1 + ^) \ - ("GmVi?) < < 00. 

(5) 

This function g{E) is continuous and smooth at E — {—Gw?/R). We define 
the entropy <S'(£') and the temperature T{E) of the system by the relations 

S{E)=\ng{E); T-\E) = f3{E) = (6) 

All the interesting thermodynamic properties of the system can be understood 
from the T{E) curve. 



6 



Consider first the case of low energies witli {—Grr? j a) < E < {—Gm^/R). 
Using ([5]) and © one can easily obtain T{E) and write it in the dimensionless 
form as 

3 r ' 



t{e) 



1 



(7) 



where we have defined t = {aT/Gm?) and e = {aE/Grn^). 

This function exhibits the peculiarities characteristic of gravitating systems. 
At the lowest energy admissible for our system, which corresponds to e = — 1 , the 
temperature t vanishes. This describes a tightly bound low temperature phase 
of the system with negligible random motion. The t{e) is clearly dominated by 
the first term of ([7]) for e ~ —1. As wc increase the energy of the system, the 
temperature increases, which is the normal behaviour for a system. This trend 
continues up to 

e = ei = -i(\/3-l)~-0.36 (8) 

at which point the t(e) curve reaches a maximum and turns around. As we 
increase the energy further the temperature decreases. The system exhibits 
negative specific heat in this range. 

Equation ([7]) is valid from the minimum energy {—Gm?/a) all the way up 
to the energy {—Gm^/R). For realistic systems, R^ a and hence this range is 
quite wide. For a small region in this range, [from {—Grn^ /a) to (— O.SGGm^/a)] 
we have positive specific heat; for the rest of the region the specific heat is 
negative. The positive specific heat region owes its existence to the nonzero 
short distance cutoff. If we set a = 0, the first term in ([7]) will vanish; we will 
have t cx (— e^^) and negative specific heat in this entire domain. 

For E > {—Gm'^/R), we have to use the second expression in ([5]) for g{E). 
In this case, we get: 



t{e) 



3[(l + e)^-f(l + fe)^] 1^-^ 
(l + e)3-(l + fe)3 



(9) 



This function, of course, matches smoothly with ^ a.t e — ~{a/ R). As we 
increase the energy, the temperature continues to decrease for a little while, 
exhibiting negative specific heat. However, this behaviour is soon halted at some 
e — £2, say. The t{e) curve reaches a minimum at this point, turns around, and 
starts increasing with increasing e. We thus enter another (high-temperature) 
phase with positive specific heat. From ^ it is clear that t ~ (l/2)e for large 
e. (Since E — {3/2)NkT for an ideal gas, we might have expected to find 
t ~ (l/3)e for our system with = 2 at high temperatures. This is indeed 
what we would have found if wc had defined our entropy as the logarithm of 
the volume of the phase space with H < E. With our definition, the energy 
of the ideal gas is actually E = [(3/2)A^ — l]kT; hence we get t — (l/2)e when 
N = 2). The form of the t(e) for (a/R) = 10"* is shown in figure [T] by the 
dashed curve. The specific heat is positive along the portions AB and CD and 
is negative along BC. 

The overall picture is now clear. Our system has two natural energy scales: 
El = (—Gm'^/a) and E2 = {—Grn^/R). For £^ > £'2, gravity is not strong 
enough to keep r < R and the system behaves like a gas confined by the con- 
tainer; we have a high temperature phase with positive specific heat. As we lower 
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Figure 1: The relation between temperature and energy for a model mimicking 
self gravitating systems. The dashed line is the result for micro-canonical en- 
semble and the solid line is for canonical ensemble. The negative specific heat 
region, BC, in the micro-canonical description is replaced by a phase transition 
in the canonical description. See text for more details 

the energy to E E2, the effects of gravity begin to be felt. For Ei < E < £^2, 

the system is unaffected by either the box or the short distance cutoff; this is 
the domain dominated entirely by gravity and we have negative specific heat. 
As we go to E El, the hard core nature of the particles begins to be felt and 
the gravity is again resisted. This gives rise to a low temperature phase with 
positive specific heat. 

We can also consider the effect of increasing R, keeping a and E fixed. Since 
we imagine the particles to be hard spheres of radius (a/2), we should only 
consider R > 2a. It is amusing to note that, if 2 < (R/a) < (a/S +1), there is 
no region of negative specific heat. As we increase R, this negative specific heat 
region appears and increasing R increases the range over which the specific heat 
is negative. Suppose a system is originally prepared with some E and R values 
such that the specific heat is positive. If we now increase R, the system may 
find itself in a region of negative specific heat. This suggests the possibility that 
an instability may be triggered in a constant energy system if its radius increases 
beyond a critical value. We will see later that this is indeed true. 
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Since systems described by canonical distribution cannot exhibit negative 
specific heat, it foUows that canonical distribution will lead to a very different 
physical picture for this range of (mean) energies Ei < E < £2- It is, therefore, 
of interest to look at our system from the point of view of canonical distribution 
by computing the partition function. In the partition function 



(10) 



the integrations over P, p and Q can be performed trivially. Omitting an overall 
constant which is unimportant, we can write the answer in the dimensionless 
form as 

^^^^ 7''" '^^'^^K^) ^''^ 

where t is the dimensionless temperature defined in ([7]) . Though this integral 
cannot be evaluated in closed form, all the limiting properties of Z{/3) can be 
easily obtained from (fTT|) . 

The integrand in (llip is large for both large and small x and reaches a 
minimum for x — Xm — (l/2i). At high temperatures, Xm < 1 and hence the 
minimum falls outside the domain of integration. The exponential contributes 
very little to the integral and we can approximate Z adequately by 



R 



R/a 



dxx 



It 



3 V a 



3a 

mi 



(12) 



On the other hand, if Xm > 1 the minimum lies between the limits of the 
integration and the exponential part of the curve dominates the integral. We 
can easily evaluate this contribution by a saddle point approach, and obtain 



r(l - 2ty^exp 



(13) 



As we lower the temperature, making Xm cross 1 from below, the contribution 
switches over from to (fT5)) . The transition is exponentially sharp. The 
critical temperature at which the transition occurs can be estimated by finding 
the temperature at which the two contributions are equal. This occurs at 



tr 



1 



1 



3 ln(i?/a) ■ 



(14) 



For t < tc, we should use p3p and for t > tc we should use (fT^ . 

Given Z((3) all thermodynamic functions can be computed. In particular, 
the mean energy of the system is given by E{(3) = —{din Z/d/S). This relation 
can be inverted to give the T{E) which can be compared with the T{E) obtained 
earlier using the micro-canonical distribution. From (fT^]) and (|13p we get. 



for t < tr and 



(15) 
(16) 
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for t > tc- Near t ~ tc, there is a rapid variation of the energy and we cannot 
use either asymptotic form. The system undergoes a phase transition ai t = tc 
absorbing a large amount of energy 

^^K^-3h^)- 

The specific heat is, of course, positive throughout the range. This is to be 
expected because canonical ensemble cannot lead to negative specific heats. 

The T — E curves obtained from the canonical (unbroken line) and micro- 
canonical (dashed line) distributions arc shown in figure[T] (For convenience, we 
have rescaled the T — E curve of the micro-canonical distribution so that e ~ 3t 
asymptotically.) At both very low and very high temperatures, the canonical 
and micro-canonical descriptions match. The crucial difference occurs at the 
intermediate energies and temperatures. Micro-canonical description predicts 
negative specific heat and a reasonably slow variation of energy with temper- 
ature. Canonical description, on the other hand, predicts a phase transition 
with rapid variation of energy with temperature. Such phase transitions are 
accompanied by large fluctuations in the energy, which is the main reason for 
the disagreement between the two descriptions [T] , [T^] , [13] . 

Numerical analysis of more realistic systems confirm all these features. Such 
systems exhibit a phase transition from the diffuse virialized phase to a core 
dominated phase when the temperature is lowered below a critical value [15j . 
The transition is very sharp and occurs at nearly constant temperature. The 
energy released by the formation of the compact core heats up the diffuse halo 
component. 



3 Isothermal sphere 

While a global maximum to the entropy does not exist in the absence of two 
cut-offs, it is, however, possible to have local extrema of entropy which are not 
global maxima. Such a configuration is described a by a Boltzmann distribu- 
tion, with the density given by p{n) oc exp[— /30(x)] where is the gravitational 
potential related to p by Poisson equation (for a formal proof see [T]). Among 
all solutions to this equation, since the one with spherical symmetry maximizes 
the entropy this configuration is usually called the isothermal sphere. The ex- 
tremum condition for the entropy, is equivalent to the differential equation for 
the gravitational potential: 

= 4^Gp,e-''['^W-^("'l (18) 

Given the solution to this equation, all other quantities can be determined. As 
we shall see, this system shows several peculiarities. 

It is convenient to introduce the length, mass and energy scale by the defi- 
nitions 

L^ = {^^Gpcfif\ Mo = 4npcLl 0^ = = ^ (19) 

where pc = p(0). All other physical variables can be expressed in terms of the 
dimensionless quantities 

x^^. nEE-^, m=^, y^/3[0-0(O)]. (20) 



10 



In terms of y{x) the isothermal equation becomes 

with the boundary condition y{Q) — y'(0) = 0. Let us consider the nature of 
solutions to this equation. 

By direct substitution, we see that n — ,to = Ix^y = 21na; satisfies 

these equations. This solution, however, is singular at the origin and hence is not 
physically admissible. The importance of this solution lies in the fact that other 
(physically admissible) solutions tend to this solution jl], [28] for large values 
of X. This asymptotic behavior of all solutions shows that the density decreases 
as (1/r^) for large r implying that the mass contained inside a sphere of radius 
r increases as M(r) oc r at large r. To find physically useful solutions, it is 
necessary to assume that the solution is cutoff at some radius R. For example, 
one may assume that the system is enclosed in a spherical box of radius R. In 
what follows, it will be assumed that the system has some cutoff radius R. 

The equation (|2ip is invariant under the transformation y y -\-a ; x —> kx 
with fc^ = e". This invariance implies that, given a solution with some value of 
y{0), we can obtain the solution with any other value of y{0) by simple rescaling. 
Therefore, only one of the two integration constants in (|2ip is really non-trivial. 
Hence it must be possible to reduce the degree of the equation from two to one 
by a judicious choice of variables [281. One such set of variables are: 

m nx^ nx'^ 
v = —; u= = . 22 

X m V 

In terms of v and w, equation IjlSp becomes 

u dv (m — 1) 



V du (u + V ~ 3) 



(23) 



The boundary conditions j/(0) = y'(0) — translate into the following: v is zero 
a.t u — 3, and {dv/du) — —5/3 at (3,0). The solution v (u) has to be obtained 
numerically: it is plotted in figure [2] as the spiraling curve. The singular points 
of this differential equation are given by the intersection of the straight lines 
u — 1 and u + V = 3 on which, the numerator and denominator of the right 
hand side of (|23p vanishes; that is, the singular point is at = 1, = 2 
corresponding to the solution n = (2/x'^), m = 2x. It is obvious from the nature 
of the equations that the solutions will spiral around the singular point. 

The nature of the solution shown in figure [5] allows us to put an interesting 
bounds on physical quantities including energy. To see this, we shall compute 
the total energy E of the isothermal sphere. The potential and kinetic energies 
are 

GM(r) dM , GM^ r» 

U = dr — — / mnxdx 

Jo r dr Lq Jq 

3M 3GM^ , , GMl3 T" 







where xq — R/Lq. The total energy is, therefore, 

E = K + U = — -2- / dx{3nx'^ - 2mnx) 
2io Jo 
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2Lo Jo dx^ ^ 



^{noxl - |mo} (25) 



where no = n{xo) and mo = m(xo). The dimensionless quantity {RE/GM'^) is 
given by 

RE 1^3, , . 

Note that the combination {RE/GM'^) is a function of{u,v) alone. Let us now 
consider the constraints on A. Suppose we specify some value for A by specifying 
R, E and M. Then such an isothermal sphere must lie on the curve 



A = 



RE 
GAP 



(27) 



which is a straight line through the point (1.5,0) with the slope A^^. On the 
other hand, since isothermal spheres must lie on the u—v curve, an isothermal 
sphere can exist only if the line in ^27\l intersects the u ~ v curve. 




Figure 2: Bound on RE /GAP for the isothermal sphere 

For large positive A (positive E) there is just one intersection. When A = 0, 
(zero energy) we still have a unique isothermal sphere. (For A = 0, equation 
([?7| is a vertical line through u = 3/2.). When A is negative (negative E), 
the line can cut the u — v curve at more than one point; thus more than one 
isothermal sphere can exist with a given value of A. [Of course, specifying 
M, i?, E individually will remove this non- uniqueness] . But as we decrease A 
(more and more negative E) the line in (|27p will slope more and more to the 
left; and when A is smaller than a critical value Ac, the intersection will cease 
to exist. Thus no isothermal sphere can exist if [RE /GA'P) is below a critical 
value Aclll This fact follows immediately from the nature oi u — v curve and 

^This derivation is due to the author I17| . It is surprising that Chandrasckhar, who has 
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equation l|27p . The value of Ac can be found from the numerical solution in 
figure. It turns out to be about (—0.335). 

The isothermal sphere has a special status as a solution to the mean field 
equations. Isothermal spheres, however, cannot exist if [RE/GM"^) < —0.335. 
Even when {RE /GA'P) > —0.335, the isothermal solution need not be stable. 
The stability of this solution can be investigated by studying the second varia- 
tion of the entropy. Such a detailed analysis shows that the following results are 
true [Hj, [H], [TTj. (i) Systems with {RE/GKP) < -0.335 cannot evolve into 
isothermal spheres. Entropy has no extremum for such systems, (ii) Systems 
with {{RE/GIvP) > -0.335) and (p(0) > 709 p(i?)) can exist in a meta-stable 
(saddle point state) isothermal sphere configuration. Here p{Q) and p{R) de- 
note the densities at the center and edge respectively. The entropy extrema 
exist but they are not local maxima, (iii) Systems with {{RE/GM"^) > —0.335) 
and {p{0) < 709 p{R)) can form isothermal spheres which are local maximum of 
entropy. 

4 An integral equation to describe nonlinear grav- 
itational clustering 

Let us next consider the gravitational clustering of a system of collision-less point 
particles in an expanding universe which poses several challenging theoretical 
questions. Though the problem can be tackled in a 'practical' manner using high 
resolution numerical simulations, such an approach hides the physical principles 
which govern the behaviour of the system. To understand the physics, it is 
necessary that we attack the problem from several directions using analytic and 
semi analytic methods. These sections will describe such attempts and will 
emphasize the semi analytic approach and outstanding issues, rather than more 
well established results. 

The expansion of the universe sets a natural length scale (called the Hubble 
radius) dn — c{d/a)^^ which is about 4000 Mpc in the current universe. In any 
region which is small compared to dfj one can set up an unambiguous coordinate 
system in which the proper coordinate of a particle r{t) = a{t)x(t) satisfies the 
Newtonian equation = — Vr^* where $ is the gravitational potential. The 
Lagrangian for such a system of particles is given by 



1 -2 

-ruirf 
2 



mimj 



(28) 



In the term 
1 



2 • 2 

a X,- 



• 2 2 

a X,- 



dx 
i—r- 

dt 



2 1 



2-2 -2 

a X,, — aax; 



dadx~ 



dt 



(29) 



we note that: (i) the total time derivative can be ignored; (ii) using d — 
— {ATrG/3)pba, the term ^frw = — (l/2)adxf = {2TTG/3)po{xf /a) can be iden- 
tified as the gravitational potential due to the uniform Friedman background of 



worked out the isothermal sphere in uv coordinates as early as 1939, missed discovering the 
energy bound shown in figure [2] Chandrasekhar 28 has the uv curve but does not over-plot 
lines of constant A. If he had done that, he would have discovered Antonov instability decades 
before Antonov did I16| . 
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density — pq/o?. Hence the Lagrangian can be expressed as 



i 

where 



ia2i2^0(t,x,) ^T-U (30) 
G rrij 2ttGpo xf 



Za ^—^ |Xj — Xj I 3 a 

is the difference between the total potential and the potential for the back ground 
Friedman universe ^frw- Varying the Lagrangian in Eq. pop with respect to 
Xi, we get the equation of motion to be: 

x + 2-x = -^V:,0 (32) 

Since is the gravitational potential generated by the perturbed mass density, 
it satisfies the equation with the source {p — ph) = pbS: 

= AT^GpbO^S (33) 

Equation (I32p and Eq. (j33p govern the nonlinear gravitational clustering in an 
expanding background. 

Usually one is interested in the evolution of the density contrast S {t, x) 
rather than in the trajectories. Since the density contrast can be expressed in 
terms of the trajectories of the particles, it should be possible to write down a 
differential equation for 5{t, x) based on the equations for the trajectories x(t) 
derived above. It is, however, somewhat easier to write down an equation for 
S]i(t) which is the spatial Fourier transform of S{t, x). To do this, we begin with 
the fact that the density p(x, t) due to a set of point particles, each of mass m, 
is given by 

p(x,i) = -^^<5z,[x-x,(<)] (34) 

where Xi(t) is the trajectory of the ith particle and Sd is the Dirac delta function. 
The density contrast S{x,t) is related to p(x,t) by 

l + ^(x,i)^^i^ = ^^5^[x-x,(<)]= fdcidD[^-Mt,cO]. (35) 
Pb N J 

In arriving at the last equality we have taken the continuum limit by: (i) re- 
placing Xi{t) by XT(i, q) where q stands for a set of parameters (like the initial 
position, velocity etc.) of a particle; for simplicity, we shall take this to be 
initial position. The subscript 'T' is just to remind ourselves that XT(i, q) is 
the trajectory of the particle, (ii) replacing {V/N) by d^q since both represent 
volume per particle. Fourier transforming both sides we get 

Skit) = J d^xc-'^''S{x,t) = J d^q exp[-ik.XT(t,q)] - i2nfSD{k) (36) 

Differentiating this expression, and using Eq. ([5^ for the trajectories give one 
can obtain an equation for (5k (see e.g. ref.[37]; Eq.lO). The structure of this 
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equation can be simplified if we use the perturbed gravitational potential (in 
Fourier space) (/)k related to 6^ by 

In terms of i/ik the exact evolution equation reads as: 

where x — XT{t,q). Of course, this equation is not 'closed'. It contains the 
velocities of the particles and their positions explicitly in the second term 
on the right and one cannot — in general — express them in simple form in 
terms of 0k- As a result, it might seem that we are in no better position than 
when we started. I will now motivate a strategy to tame this term in order to 
close this equation. This strategy depends on two features: 

(1) First, extremely nonlinear structures do not contribute to the right hand 
side of Eq. ([S5|) though, of course, they contribute individually to the two terms. 
More precisely, the right hand side of Eq. ([38|) will lead to a density contrast 
that will scale as i5k cx; k'^ if originally — in linear theory — 6k oc k" with n > 2 
as fc 0. This leads to a P oc fc^ tail in the power spectrum [551 [57]. (We will 
provide a derivation of this result in the next section.) 

(ii) Second, we can use Zeldovich approximation to evaluate this term, once 
the above fact is realized. It is well known that, when the density contrasts are 
small, it grows as 5k oc a in the linear limit. One can easily show that such a 
growth corresponds to particle displacements of the form 

XT(a,q) =q-aVi/^(q); = (4^Gpo)~V (39) 

A useful approximation to describe the quasi linear stages of clustering is ob- 
tained by using the trajectory in Eq. (j39p as an ansatz valid even at quasi linear 
epochs. In this approximation, (called Zeldovich approximation), the velocities 
x can be expressed in terms of the initial gravitational potential. 

We now combine the two results mentioned above to obtain a closure con- 
dition for our dynamical equation. At any given moment of time we can divide 
the particles in the system into three different sets. First, there are particles 
which are already a part of virialized cluster in the non linear regime. Second 
set is made of particles which are completely unbound and are essentially con- 
tributing to power spectrum at the linear scales. The third set is made of all 
particles which cannot be put into either of these two baskets. Of these three, 
we know that the first two sets of particles do not contribute significantly to 
the right hand side of Eq. psp so we will not incur any serious error in ignoring 
these particles in computing the right hand side. For the description of particles 
in the third set, the Zeldovich approximation should be fairly good. In fact, we 
can do slightly better than the standard Zeldovich approximation. Wc note that 
in Eq. ([39]) the velocities were taken to be proportional to the gradient of the 
initial gravitational potential. We can improve on this ansatz by taking the ve- 
locities to be given by the gradient of the instantaneous gravitational potential 



k + 4-</)k 
a 



1 f d^p 

(2^ 



2a2 
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which has the effect of incorporating the influence of particles in bound clusters 
on the rest of the particles to certain extent. 

Given this ansatz, it is straightforward to obtain a closed integro-differential 
equation for the gravitational potential. Direct calculation shows that the grav- 
itational potential is described by the closed integral equation: (The details of 
this derivation can be found in ref. [27j and will not be repeated here.) 



a 6a'' 



(2^ 



'ik+pfpik-p 



7;P - 5 



k p 



(40) 



This equation provides a powerful method for analyzing non linear clustering 
since estimating Eq. (|38p by Zeldovich approximation has a very large domain 
of applicability. In the next two sections, I will use this equation to study the 
transfer of power in gravitational clustering. 



5 Inverse cascade in non linear gravitational clus- 
tering: The /c^ tail 

There is an interesting and curious result which is characteristic of gravitational 
clustering that can be obtained directly from our Eq. (|40|) . Consider an initial 
power spectrum which has very little power at large scales; more precisely, we 
shall assume that P(k) oc k" with n > 4 for small k (i.e, the power dies faster 
than k'^ for small k) . If these large spatial scales are described by linear theory — 
as one would have normally expected — then the power at these scales can only 
grow as P oc a^fc" and it will always be sub dominant to fc^. It turns out that 
this conclusion is incorrect. As the system evolves, small scale nonlinearities will 
develop in the system and — if the large scales have too little power intrinsically 
— then the long wavelength power will soon be dominated by the "fc'*-tail" of 
the short wavelength power arising from the nonlinear clustering. This is a 
purely non linear effect which we shall now describe. (This result is known in 
literature [22l [27] but we will derive it from the formalism developed in the last 
section which adds fresh insight.) 

A formal way of obtaining the fc^ tail is to solve Eq. for long wavelengths; 
i.e. near k = 0. Writing (p^ = 4''^^ + 4>k^ + where (f)'^'' — cj)^'' is the time 

(2) 

independent gravitational potential in the linear theory and cjy^ is the next 
order correction, we get from Eq. (HO]), the equation 

where 5(k,p) = [(7/8)A:2 + (3/2)p2_5(k-p/fc))2. The solution to this equation is 
the sum of a solution to the homogeneous part [which decays as </) oc a~'^ oc t~^^^ 
giving (j) oc t^^/"^] and a particular solution which grows as a. Ignoring the 

(2) 

decaying mode at late times and taking (jr^ — aC\^ one can determine Ck from 
the above equation. Plugging it back, we find the lowest order correction to be. 
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Figure 3: The transfer of power to long wavelengths forming a tail is illus- 
trated using simulation results. Power is injected in the form of a narrow peak 
at L = 8. Note that the y— axis is (A/a) so that there will be no change of 
shape of the power spectrum under linear evolution with A cx a. As time goes 
on a /c** tail is generated purely due to nonlinear coupling between the modes. 
(Figure adapted from ref.|30j.) 



Near k ~ 0, we have 



^(2) 



2a 
2LHf 



(2^ 



7 3 5(k.p)^ 
s'' +2^ fc2 



1267r2i72 



L)|2 



(43) 



which is independent of k to the lowest order. Correspondingly the power 
spectrum for density Ps{k) cx a?k'^P^{k) cx a^k'^ in this order. 

The generation of long wavelength fc^ tail is easily seen in simulations if one 
starts with a power spectrum that is sharply peaked in |k|. Figure [3] (adapted 
from [30]) shows the results of such a simulation. The y-axis is [A(fc)/a(i)] 
where A2(fc) = k^Pjl-Tx'^ is the power per logarithmic band in k. In linear 
theory A cx a and this quantity should not change. The curves labelled by 
a = 0.12 to a = 20.0 show the effects of nonlinear evolution, especially the 
development of fc^ tail. (Actually one can do better. The formation can also 
reproduce the sub-harmonic at i ~ 4 seen in Fig. [3] and other details; see ref. 

m) 
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6 Analogue of Kolmogorov spectrum for gravi- 
tational clustering 

If power is injected at some scale L into an ordinary viscous fluid, it cascades 
down to smaller scales because of the non linear coupling between different 
modes. The resulting power spectrum, for a wide range of scales, is well ap- 
proximated by the Kolmogorov spectrum which plays a key, useful, role in the 
study of fluid turbulence. It is possible to obtain the form of this spectrum 
from fairly simple and general considerations though the actual equations of 
fluid turbulence are intractably complicated. Let us now consider the corre- 
sponding question for non linear gravitational clustering. If power is injected at 
a given length scale very early on, how does the dynamical evolution transfers 
power to other scales at late times? In particular, does the non linear evolution 
lead to an analogue of Kolmogorov spectrum with some level of universality, in 
the case of gravitational interactions? 

Surprisingly, the answer is "yes" , even though normal fluids and collisionlcss 
self gravitating particles constitute very different physical systems. If power is 
injected at a given scale L = 27r/fco then the gravitational clustering transfers 
the power to both larger and smaller spatial scales. At large spatial scales the 
power spectrum goes as P{k) cx fc^ as soon as non linear coupling becomes 
important. We have already seen this result in the previous section. More 
interestingly, the cascading of power to smaller scales leads to a universal pattern 
at late times just as in the case of fluid turbulence. This is because, Eq. ((40|) 
admits solutions for the gravitational potential of the form 0k(^) = 
at late times when the initial condition is irrelevant; here F{t) satisfies a non 
linear differential equation and I?(k) satisfies an integral equation. One can 
analyze the relevant equations analytically as well as verified the conclusions by 
numerical simulations. This study (the details of which can be found in ref.|31j) 
confirms that non linear gravitational clustering does lead to a universal power 
spectrum at late times if the power is injected at a given scale initially. (In 
cosmology there is very little motivation to study the transfer of power by itself 
and most of the numerical simulations in the past concentrated on evolving 
broad band initial power spectrum. So this result was missed out.) 

Our aim is to look for late time scale free evolution of the system exploiting 
the fact Eq. ipp)) allows self similar solutions of the form 4)-k_{t) = F{t)D{'k). 
Substituting this ansatz into Eq. (|40p we obtain two separate equations for 
F{t) and -D(k). It is also convenient at this stage to use the expansion factor 
a{t) — {t/tf)Y/^ of the matter dominated universe as the independent variable 
rather than the cosmic time t. Then simple algebra shows that the governing 
equations are 

a—+'L—~-F^ (44) 
da"^ 2 da 

and 

^IJ (^^^k+p^ik-pe(k, P) (45) 

Equation (|44p governs the time evolution while Eq. governs the shape of 
the power spectrum. (The separation ansatz, of course, has the scaling freedom 
F — > nF,D {l/fi)D which will change the right hand side of Eq. to 
—fjiF^ and the left hand side of Eq. (|^5)) to jiH^D]^. But, as to be expected, our 
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Figure 4: The solution to Eq. ((46|) plotted in g — a plane. The function g{a) 
asymptotically approaches unity with oscillations which are represented by the 
spiral in the right panel. The different curves in the left panel corresponds to 
the rescaling freedom in the initial conditions. One fiducial curve which was 
used to model the simulation is shown in the red. For more details, see ref. [SI] . 

results will be independent of /i; so we have set it to unity). Our interest lies 
in analyzing the solutions of Eq. (j44p subject to the initial conditions F{ai) — 
constant, {dF/da)i at some small enough a — ai. 

Inspection shows that Eq. ([33]) has the exact solution F{a) — (3/2)a^^. 
This, of course, is a special solution and will not satisfy the relevant initial con- 
ditions. However, Eq. ([Uj) fortunately belongs to a class of non linear equations 
which can be mapped to a homologous system. In such cases, the special power 
law solutions will arise as the asymptotic limit. (The example well known to 
astronomers is that of isothermal sphere [2H]- Our analysis below has a close 
parallel.) To find the general behaviour of the solutions to Eq. ([Hj), we will make 
the substitution F{a) ~ (3/2)a^^g(a) and change the independent variable from 
aio q — logo. Then Eq. (I3i|) reduces to the form 



This represents a particle moving in a potential V{g) ~ {l/2)g'^ — (3/4)g^ under 
friction. For our initial conditions the motion will lead the "particle" to asymp- 
totically come to rest at the stable minimum at 5 = 1 with damped oscillations. 
In other words, F{a) —>■ (3/2)a^^ for large a showing this is indeed the asymp- 
totic solution. From the Poisson equation, it follows that fc^(/>k oc ((5k/a) so that 
(5k(a) oc g{a)k^D(h) giving a direct physical meaning to the function g{a) as 




(46) 
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° 10 100 =■ 10 100 

L=2 Tr/k L-2 n/k 

Figure 5: Left panel: The results of the numerical simulation with an initial 
power spectrum which is a Gaussian peaked at L — 24. The y-axis gives 
where A| = k^P/2T:'^ is the power per logarithmic band. The evolution gen- 
erates a well known fc^ tail at large scales (see for example, [30]) and leads 
to cascading of power to small scales. Right panel: The simulation data is 
re-expressed by factoring out the time evolution function g{a) obtained by 
integrating Eq. The fact that the curves fall nearly on top of each 

other shows that the late time evolution is scale free and described by the 
ansatz discussed in the text. The rescaled spectrum is very well described by 
P{k) cx fc^° '*/(l + (k/ko)'^'^) which is shown by the, completely overlapping, 
broken blue curve. For more details, see ref. [31j. 
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the growth factor for the density contrast. The asymptotic Umit {g ~ 1) cor- 
responds to to a rather trivial case of 6^ becoming independent of time. What 
will be more interesting — and accessible in simulations — will be the approach 
to this asymptotic solution. To obtain this, we introduce the variable 

u = g + 2(^) (47) 



dq 

so that our system becomes homologous. It can be easily shown that we now 
get the first order form of the autonomous system to be 

du 6q(q — 1) , 
— = — — 48) 

dg u- g 

The critical points of the system are at (0, 0) and (1,1). Standard analysis shows 
that: (i) the (0,0) is an unstable critical point and the second one (1, 1) is the 
stable critical point; (ii) for our initial conditions the solution spirals around the 
stable critical point. 

Figure [3] (from ref. [3T]) describes the solution in the g — a plane. The 
g{a) curves clearly approach the asymptotic value of g « 1 with superposed 
oscillations. The different curves in Fig. ^ are for different initial values which 
arise from the scaling freedom mentioned earlier. (The thick red line correspond 
to the initial conditions used in the simulations described below.). The solu- 
tion g{a) describes the time evolution and solves the problem of determining 
asymptotic time evolution. 

To test the correctness of these conclusions, we performed a high resolution 
simulation using the TreePM method [32l [33] and its parallel version [34] with 
128^ particles on a 128'^ grid. Details about the code parameters can be found 
in [33] . The initial power spectrum P{k) was chosen to be a Gaussian peaked at 
the scale of kp — 2tt/ Lp with Lp — 24 grid lengths and with a standard deviation 
Ak = 2tt/ Lbox, where Lbox — 128 is the size of one side of the simulation volume. 
The amplitude of the peak was taken such that A/i„ [kp — 27r/Lp, a — 0.25) — 1. 

The late time evolution of the power spectrum (in terms of = k^P{k) /27r^ 
where P = \Sk\'^ is the power spectrum of density fluctuations) obtained from the 
simulations is shown in FigIS] (left panel) . In the right panel, we have rescaled 
the Afe, using the appropriate solution g{a). The fact that the curves fall on 
top of each other shows that the late time evolution indeed sales as g(a) within 
numerical accuracy. A reasonably accurate fit for g(a) at late times used in this 
figure is given by g{a) oc a(l — 0.3 In a). The key point to note is that the asymp- 
totic time evolution is essentially S{a) oc a except for a logarithmic correction, 
even in highly nonlinear scales. (This was first noticed from somewhat lower 
resolution simulations in [3D].). Since the evolution at linear scales is always 
(5 oc a, this allows for a form invariant evolution of power spectrum at all scales. 
Gravitational clustering evolves towards this asymptotic state. 

To the lowest order of accuracy, the power spectrum at this range of scales 
is approximated by the mean index n w — 1 with P{k) oc k~^. A better fit to 
the power spectrum in FigEJis given by 

l + (k/koy-^ ko 

This fit is shown by the broken blue line in the figure which completely overlaps 
with the data and is barely visible. (Note that this fit is applicable only at 
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L < Lp since the fc"* tail will dominate scales to the right of the initial peak; 
see the discussion in [ID])- At nonlinear scales P{k) oc making flat, as 
seen in Fig. \E\ (This is not a numerical artifact and we have sufficient dynamic 
range in the simulation to ascertain this.) At quasi linear scales P{k) oc k^'^-'^. 
The effective index of the power spectrum varies between —3 and —0.4 in this 
range of scales. 

What could be a possible interpretation for this behaviour ? It is difhcult 
to provide a simple but precise answer but one possible line of reasoning is as 
follows: In the case of viscous fluids, the energy is dissipated at the smallest 
scales as heat. In steady state, energy cannot accumulate at any intermediate 
scale and hence the rate of flow of energy from one scale to the next (lower) 
scale must be a constant. This constancy immediately leads to Kolmogorov 
spectrum. In the case of gravitating particles, there is no dissipation and each 
scale will evolve towards virial equilibrium. At any given time t, the power 
would have cascaded down only up to some scale Iminit) which it self, of course, 
is a decreasing function of time. So, at time t we expect very little power for 
1 < klmin{t) and a k^ tail for kLp < 1, say. The really interesting band is 
between Imin and Lp. 

To understand this band, let us recall that the Lagrangian in Eq. ((50)) leads to 
the time dependent Hamiltonian is H(p, x, t) = ^[p^/2ma^ + C/]. The evolution 
of the energy in the system is governed by the equation dH/da = {dH/da)p^x. 
It is clear from Eq.^ that {dU/da)p^^ = -U/a while (9T/(9a)p,x = ~2T/a. 
Hence the time evolution of the total energy H = E oi the system is described 

by 

dE 1, , 2E U E T 

— = {2T + U)^ = (50) 

da a a a a a 

In the continuum limit, ignoring the infinite self-energy term, the potential 
energy can be written as: 

2a ./ ./ x-y ^ ^ 



d{a^E) _ ^^^^_Gpla f ^ f ^ S{K,a)S{y,a) 



Hence 

'^'"^^^ = -a'U = ^ / d^x / rf3 rvrrpiA:^ (52) 
da 2 7 7 |x-y| ^ ^ 

The ensemble average of the right hand side, per unit proper volume will be 

2Va^ J |x - y| J a^k^ Jq k 

where V is the comoving volume. 

When a particular scale is virialized, we expect £ ~ constant at that scale in 
comoving coordinates. That is, we would expect the energy density in Eq. (|53p 
would have reached equipartition and contribute same amount per logarithmic 
band of scales in the intermediate scales between Imm and Lpeak- This requires 
P{k) cx a^/fc which is essentially what we found from simulations. The time de- 
pendence of P is essentially P oc a except for a logarithmic correction. Similarly 
the scale dependence is P oc fc^^ which is indeed a good fit to the simulation 
results. The flattening of the power at small scales, modeled by the more precise 
fitting function in Eq. ()49p . can be understood from the fact that, equipartition 
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is not yet achieved at smaller scales. The same result holds for kinetic energy 
if the motion is dominated by scale invariant radial flows [301 135j . Our result 
suggests that gravitational power transfer evolves towards this equipartition. 
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